!*************************************************************************
!Copyright (c) 1980-2007, Jorge Nocedal.
!
!Contributors:
!    * Sergey Bochkanov (ALGLIB project). Translation from FORTRAN to
!      pseudocode.
!
!This software is freely available for educational or commercial  purposes.
!We expect that all publications describing work using this software  quote
!at least one of the references given below:
!    * J. Nocedal. Updating  Quasi-Newton  Matrices  with  Limited  Storage
!      (1980), Mathematics of Computation 35, pp. 773-782.
!    * D.C. Liu and J. Nocedal. On the  Limited  Memory  Method  for  Large
!      Scale  Optimization  (1989),  Mathematical  Programming  B,  45,  3,
!      pp. 503-528.
!*************************************************************************)
module lbfgs

!
!This routines must be defined by you:
!subroutine FuncGrad(X : TReal1DArray; var F : Double; var G : TReal1DArray);
!

contains

!(*************************************************************************
!        LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
!                          JORGE NOCEDAL

!The subroutine minimizes function F(x) of N arguments by  using  a  quasi-
!Newton method (LBFGS scheme) which is optimized to use  a  minimum  amount
!of memory.

!The subroutine generates the approximation of an inverse Hessian matrix by
!using information about the last M steps of the algorithm  (instead of N).
!It lessens a required amount of memory from a value  of  order  N^2  to  a
!value of order 2*N*M.

!This subroutine uses the FuncGrad subroutine which calculates the value of
!the function F and gradient G in point X. The programmer should define the
!FuncGrad subroutine by himself.  It  should  be  noted that the subroutine
!doesn't need to waste time for memory allocation of array G,  because  the
!memory is allocated in calling the subroutine. Setting a dimension of array
!G  each  time  when  calling  a  subroutine  will excessively slow down an
!algorithm.

!The programmer could also redefine the LBFGSNewIteration subroutine  which
!is called on each new step. The current point X, the function value F  and
!the  gradient  G  are  passed  into  this  subroutine. It is reasonable to
!redefine the subroutine for better debugging, for  example,  to  visualize
!the solution process.

!Input parameters:
!    N   -   problem dimension. N>0
!    M   -   number of corrections in the BFGS scheme of Hessian
!            approximation update. Recommended value:  3<=M<=7. The smaller
!            value causes worse convergence, the bigger will  not  cause  a
!            considerably better convergence, but will cause a fall in  the
!            performance. M<=N.
!    X   -   initial solution approximation.
!            Array whose index ranges from 1 to N.
!    EpsG -  positive number which  defines  a  precision  of  search.  The
!            subroutine finishes its work if the condition ||G|| < EpsG  is
!            satisfied, where ||.|| means Euclidian norm, G - gradient, X -
!            current approximation.
!    EpsF -  positive number which  defines  a  precision  of  search.  The
!            subroutine finishes its work if on iteration  number  k+1  the
!            condition |F(k+1)-F(k)| <= EpsF*max{|F(k)|, |F(k+1)|, 1}    is
!            satisfied.
!    EpsX -  positive number which  defines  a  precision  of  search.  The
!            subroutine finishes its work if on iteration number k+1    the
!            condition |X(k+1)-X(k)| <= EpsX is fulfilled.
!    MaxIts- maximum number of iterations. If MaxIts=0, the number of
!            iterations is unlimited.
!
!Output parameters:
!    X   -   solution approximation. Array whose index ranges from 1 to N.
!    Info-   a return code:
!                    * -1 wrong parameters were specified,
!                    * 0 interrupted by user,
!                    * 1 relative function decreasing is less or equal to EpsF,
!                    * 2 step is less or equal EpsX,
!                    * 4 gradient norm is less or equal to EpsG,
!                    * 5 number of iterations exceeds MaxIts.
!
!FuncGrad routine description. User-defined.
!Input parameters:
!    X   -   array whose index ranges from 1 to N.
!Output parameters:
!    G   -   function gradient.
!            Array whose index ranges from 1 to N.
!The memory for array G has already been allocated in the calling subroutine,
!and it isn't necessary to allocate it in the FuncGrad subroutine.
!*************************************************************************)

subroutine LBFGSMinimize(N, M, X, EpsG, EpsF, EpsX, MaxIts, Info)
   implicit none
   integer, intent(in) :: N, M, MaxIts

     const N : Integer;
     const M : Integer;
     var X : TReal1DArray;
     const EpsG : Double;
     const EpsF : Double;
     const EpsX : Double;
     const MaxIts : Integer;
     var Info : Integer);
var
    W : TReal1DArray;
    F : Double;
    FOld : Double;
    TF : Double;
    TXNorm : Double;
    V : Double;
    XOld : TReal1DArray;
    TX : TReal1DArray;
    G : TReal1DArray;
    Diag : TReal1DArray;
    TA : TReal1DArray;
    FINISH : Boolean;
    GNORM : Double;
    STP1 : Double;
    FTOL : Double;
    STP : Double;
    YS : Double;
    YY : Double;
    SQ : Double;
    YR : Double;
    BETA : Double;
    XNORM : Double;
    ITER : Integer;
    NFUN : Integer;
    POINT : Integer;
    ISPT : Integer;
    IYPT : Integer;
    MAXFEV : Integer;
    BOUND : Integer;
    NPT : Integer;
    CP : Integer;
    I : Integer;
    NFEV : Integer;
    INMC : Integer;
    IYCN : Integer;
    ISCN : Integer;
    XTOL : Double;
    GTol : Double;
    STPMIN : Double;
    STPMAX : Double;
begin
    SetLength(W, N*(2*M+1)+2*M+1);
    SetLength(G, N+1);
    SetLength(XOld, N+1);
    SetLength(TX, N+1);
    SetLength(Diag, N+1);
    SetLength(TA, N+1);
    FuncGrad(X, F, G);
    FOld := F;
    ITER := 0;
    Info := 0;
    if (N<=0) or (M<=0) or (M>N) or (EpsG<0) or (EpsF<0) or (EpsX<0) or (MaxIts<0) then
    begin
        Info := -1;
        Exit;
    end;
    NFUN := 1;
    POINT := 0;
    FINISH := False;
    I:=1;
    while I<=N do
    begin
        DIAG[I] := 1;
        Inc(I);
    end;
    XTOL := 100*MachineEpsilon;
    GTol := 0.9;
    STPMIN := Power(10, -20);
    STPMAX := Power(10, 20);
    ISPT := N+2*M;
    IYPT := ISPT+N*M;
    I:=1;
    while I<=N do
    begin
        W[ISPT+I] := -G[I]*DIAG[I];
        Inc(I);
    end;
    GNORM := Sqrt(LBFGSDotProduct(N, G, 1, G, 1));
    if GNORM<=EpsG then
    begin
        Info := 4;
        Exit;
    end;
    STP1 := 1/GNORM;
    FTOL := 0.0001;
    MAXFEV := 20;
    while True do
    begin
        APVMove(@XOld[0], 1, N, @X[0], 1, N);
        ITER := ITER+1;
        INFO := 0;
        BOUND := ITER-1;
        if ITER<>1 then
        begin
            if ITER>M then
            begin
                BOUND := M;
            end;
            YS := LBFGSDotProduct(N, W, IYPT+NPT+1, W, ISPT+NPT+1);
            YY := LBFGSDotProduct(N, W, IYPT+NPT+1, W, IYPT+NPT+1);
            I:=1;
            while I<=N do
            begin
                DIAG[I] := YS/YY;
                Inc(I);
            end;
            CP := POINT;
            if POINT=0 then
            begin
                CP := M;
            end;
            W[N+CP] := 1/YS;
            I:=1;
            while I<=N do
            begin
                W[I] := -G[I];
                Inc(I);
            end;
            CP := POINT;
            I:=1;
            while I<=BOUND do
            begin
                CP := CP-1;
                if CP=-1 then
                begin
                    CP := M-1;
                end;
                SQ := LBFGSDotProduct(N, W, ISPT+CP*N+1, W, 1);
                INMC := N+M+CP+1;
                IYCN := IYPT+CP*N;
                W[INMC] := W[N+CP+1]*SQ;
                LBFGSLinComb(N, -W[INMC], W, IYCN+1, W, 1);
                Inc(I);
            end;
            I:=1;
            while I<=N do
            begin
                W[I] := DIAG[I]*W[I];
                Inc(I);
            end;
            I:=1;
            while I<=BOUND do
            begin
                YR := LBFGSDotProduct(N, W, IYPT+CP*N+1, W, 1);
                BETA := W[N+CP+1]*YR;
                INMC := N+M+CP+1;
                BETA := W[INMC]-BETA;
                ISCN := ISPT+CP*N;
                LBFGSLinComb(N, BETA, W, ISCN+1, W, 1);
                CP := CP+1;
                if CP=M then
                begin
                    CP := 0;
                end;
                Inc(I);
            end;
            I:=1;
            while I<=N do
            begin
                W[ISPT+POINT*N+I] := W[I];
                Inc(I);
            end;
        end;
        NFEV := 0;
        STP := 1;
        if ITER=1 then
        begin
            STP := STP1;
        end;
        I:=1;
        while I<=N do
        begin
            W[I] := G[I];
            Inc(I);
        end;
        LBFGSMCSRCH(N, X, F, G, W, ISPT+POINT*N+1, STP, FTOL, XTOL, MAXFEV, INFO, NFEV, DIAG, GTOL, STPMIN, STPMAX);
        if INFO<>1 then
        begin
            if Info=0 then
            begin
                Info := -1;
                Exit;
            end;
        end;
        NFUN := NFUN+NFEV;
        NPT := POINT*N;
        I:=1;
        while I<=N do
        begin
            W[ISPT+NPT+I] := STP*W[ISPT+NPT+I];
            W[IYPT+NPT+I] := G[I]-W[I];
            Inc(I);
        end;
        POINT := POINT+1;
        if POINT=M then
        begin
            POINT := 0;
        end;
        if (ITER>MaxIts) and (MaxIts>0) then
        begin
            Info := 5;
            Exit;
        end;
        LBFGSNewIteration(X, F, G);
        GNORM := Sqrt(LBFGSDotProduct(N, G, 1, G, 1));
        if GNORM<=EpsG then
        begin
            Info := 4;
            Exit;
        end;
        TF := Max(AbsReal(FOld), Max(AbsReal(F), 1.0));
        if FOld-F<=EpsF*TF then
        begin
            Info := 1;
            Exit;
        end;
        APVMove(@TX[0], 1, N, @XOld[0], 1, N);
        APVSub(@TX[0], 1, N, @X[0], 1, N);
        XNORM := Sqrt(LBFGSDotProduct(N, X, 1, X, 1));
        TXNorm := Max(XNORM, Sqrt(LBFGSDotProduct(N, XOld, 1, XOld, 1)));
        TXNorm := Max(TXNorm, 1.0);
        V := Sqrt(LBFGSDotProduct(N, TX, 1, TX, 1));
        if V<=EpsX then
        begin
            Info := 2;
            Exit;
        end;
        FOld := F;
        APVMove(@XOld[0], 1, N, @X[0], 1, N);
    end;
end;


subroutine LBFGSLinComb(const N : Integer;
     const DA : Double;
     const DX : TReal1DArray;
     SX : Integer;
     var DY : TReal1DArray;
     SY : Integer);
var
    fx : Integer;
    fy : Integer;
begin
    fx := sx+n-1;
    fy := sy+n-1;
    APVAdd(@dy[0], sy, fy, @dx[0], sx, fx, da);
end;


function LBFGSDotProduct(const n : Integer;
     const dx : TReal1DArray;
     sx : Integer;
     const dy : TReal1DArray;
     sy : Integer):Double;
var
    v : Double;
    fx : Integer;
    fy : Integer;
begin
    fx := sx+n-1;
    fy := sy+n-1;
    v := APVDotProduct(@dx[0], sx, fx, @dy[0], sy, fy);
    Result := v;
end;


subroutine LBFGSMCSRCH(const N : Integer;
     var X : TReal1DArray;
     var F : Double;
     var G : TReal1DArray;
     const S : TReal1DArray;
     SStart : Integer;
     var STP : Double;
     const FTOL : Double;
     const XTOL : Double;
     const MAXFEV : Integer;
     var INFO : Integer;
     var NFEV : Integer;
     var WA : TReal1DArray;
     const GTOL : Double;
     const STPMIN : Double;
     const STPMAX : Double);
var
    INFOC : Integer;
    J : Integer;
    BRACKT : Boolean;
    STAGE1 : Boolean;
    DG : Double;
    DGM : Double;
    DGINIT : Double;
    DGTEST : Double;
    DGX : Double;
    DGXM : Double;
    DGY : Double;
    DGYM : Double;
    FINIT : Double;
    FTEST1 : Double;
    FM : Double;
    FX : Double;
    FXM : Double;
    FY : Double;
    FYM : Double;
    P5 : Double;
    P66 : Double;
    STX : Double;
    STY : Double;
    STMIN : Double;
    STMAX : Double;
    WIDTH : Double;
    WIDTH1 : Double;
    XTRAPF : Double;
    ZERO : Double;
    MyTemp : Double;
begin
    SStart := SStart-1;
    P5 := 0.5;
    P66 := 0.66;
    XTRAPF := 4.0;
    ZERO := 0;
    FuncGrad(X, F, G);
    INFOC := 1;
    INFO := 0;
    if (N<=0) or (STP<=0) or (FTOL<0) or (GTOL<ZERO) or (XTOL<ZERO) or (STPMIN<ZERO) or (STPMAX<STPMIN) or (MAXFEV<=0) then
    begin
        Exit;
    end;
    DGINIT := 0;
    J:=1;
    while J<=N do
    begin
        DGINIT := DGINIT+G[J]*S[J+SStart];
        Inc(J);
    end;
    if DGINIT>=0 then
    begin
        Exit;
    end;
    BRACKT := False;
    STAGE1 := True;
    NFEV := 0;
    FINIT := F;
    DGTEST := FTOL*DGINIT;
    WIDTH := STPMAX-STPMIN;
    WIDTH1 := WIDTH/P5;
    J:=1;
    while J<=N do
    begin
        WA[J] := X[J];
        Inc(J);
    end;
    STX := 0;
    FX := FINIT;
    DGX := DGINIT;
    STY := 0;
    FY := FINIT;
    DGY := DGINIT;
    while True do
    begin
        if BRACKT then
        begin
            if STX<STY then
            begin
                STMIN := STX;
                STMAX := STY;
            end
            else
            begin
                STMIN := STY;
                STMAX := STX;
            end;
        end
        else
        begin
            STMIN := STX;
            STMAX := STP+XTRAPF*(STP-STX);
        end;
        if STP>STPMAX then
        begin
            STP := STPMAX;
        end;
        if STP<STPMIN then
        begin
            STP := STPMIN;
        end;
        if BRACKT and ((STP<=STMIN) or (STP>=STMAX)) or (NFEV>=MAXFEV-1) or (INFOC=0) or BRACKT and (STMAX-STMIN<=XTOL*STMAX) then
        begin
            STP := STX;
        end;
        J:=1;
        while J<=N do
        begin
            X[J] := WA[J]+STP*S[J+SStart];
            Inc(J);
        end;
        FuncGrad(X, F, G);
        INFO := 0;
        NFEV := NFEV+1;
        DG := 0;
        J:=1;
        while J<=N do
        begin
            DG := DG+G[J]*S[J+SStart];
            Inc(J);
        end;
        FTEST1 := FINIT+STP*DGTEST;
        if BRACKT and ((STP<=STMIN) or (STP>=STMAX)) or (INFOC=0) then
        begin
            INFO := 6;
        end;
        if (STP=STPMAX) and (F<=FTEST1) and (DG<=DGTEST) then
        begin
            INFO := 5;
        end;
        if (STP=STPMIN) and ((F>FTEST1) or (DG>=DGTEST)) then
        begin
            INFO := 4;
        end;
        if NFEV>=MAXFEV then
        begin
            INFO := 3;
        end;
        if BRACKT and (STMAX-STMIN<=XTOL*STMAX) then
        begin
            INFO := 2;
        end;
        if (F<=FTEST1) and (AbsReal(DG)<=-GTOL*DGINIT) then
        begin
            INFO := 1;
        end;
        if INFO<>0 then
        begin
            Exit;
        end;
        MyTemp := FTOL;
        if GTOL<FTOL then
        begin
            MyTemp := GTOL;
        end;
        if STAGE1 and (F<=FTEST1) and (DG>=MyTemp*DGINIT) then
        begin
            STAGE1 := False;
        end;
        if STAGE1 and (F<=FX) and (F>FTEST1) then
        begin
            FM := F-STP*DGTEST;
            FXM := FX-STX*DGTEST;
            FYM := FY-STY*DGTEST;
            DGM := DG-DGTEST;
            DGXM := DGX-DGTEST;
            DGYM := DGY-DGTEST;
            LBFGSMCSTEP(STX, FXM, DGXM, STY, FYM, DGYM, STP, FM, DGM, BRACKT, STMIN, STMAX, INFOC);
            FX := FXM+STX*DGTEST;
            FY := FYM+STY*DGTEST;
            DGX := DGXM+DGTEST;
            DGY := DGYM+DGTEST;
        end
        else
        begin
            LBFGSMCSTEP(STX, FX, DGX, STY, FY, DGY, STP, F, DG, BRACKT, STMIN, STMAX, INFOC);
        end;
        if BRACKT then
        begin
            if AbsReal(STY-STX)>=P66*WIDTH1 then
            begin
                STP := STX+P5*(STY-STX);
            end;
            WIDTH1 := WIDTH;
            WIDTH := AbsReal(STY-STX);
        end;
    end;
end;


subroutine LBFGSMCSTEP(var STX : Double;
     var FX : Double;
     var DX : Double;
     var STY : Double;
     var FY : Double;
     var DY : Double;
     var STP : Double;
     const FP : Double;
     const DP : Double;
     var BRACKT : Boolean;
     const STMIN : Double;
     const STMAX : Double;
     var INFO : Integer);
var
    BOUND : Boolean;
    GAMMA : Double;
    P : Double;
    Q : Double;
    R : Double;
    S : Double;
    SGND : Double;
    STPC : Double;
    STPF : Double;
    STPQ : Double;
    THETA : Double;
begin
    INFO := 0;
    if BRACKT and ((STP<=Min(STX, STY)) or (STP>=Max(STX, STY))) or (DX*(STP-STX)>=0) or (STMAX<STMIN) then
    begin
        Exit;
    end;
    SGND := DP*(DX/AbsReal(DX));
    if FP>FX then
    begin
        INFO := 1;
        BOUND := True;
        THETA := 3*(FX-FP)/(STP-STX)+DX+DP;
        S := Max(AbsReal(THETA), Max(AbsReal(DX), AbsReal(DP)));
        GAMMA := S*Sqrt(Sqr(THETA/S)-DX/S*(DP/S));
        if STP<STX then
        begin
            GAMMA := -GAMMA;
        end;
        P := GAMMA-DX+THETA;
        Q := GAMMA-DX+GAMMA+DP;
        R := P/Q;
        STPC := STX+R*(STP-STX);
        STPQ := STX+DX/((FX-FP)/(STP-STX)+DX)/2*(STP-STX);
        if AbsReal(STPC-STX)<AbsReal(STPQ-STX) then
        begin
            STPF := STPC;
        end
        else
        begin
            STPF := STPC+(STPQ-STPC)/2;
        end;
        BRACKT := True;
    end
    else
    begin
        if SGND<0 then
        begin
            INFO := 2;
            BOUND := False;
            THETA := 3*(FX-FP)/(STP-STX)+DX+DP;
            S := Max(AbsReal(THETA), Max(AbsReal(DX), AbsReal(DP)));
            GAMMA := S*SQRT(Sqr(THETA/S)-DX/S*(DP/S));
            if STP>STX then
            begin
                GAMMA := -GAMMA;
            end;
            P := GAMMA-DP+THETA;
            Q := GAMMA-DP+GAMMA+DX;
            R := P/Q;
            STPC := STP+R*(STX-STP);
            STPQ := STP+DP/(DP-DX)*(STX-STP);
            if AbsReal(STPC-STP)>AbsReal(STPQ-STP) then
            begin
                STPF := STPC;
            end
            else
            begin
                STPF := STPQ;
            end;
            BRACKT := True;
        end
        else
        begin
            if AbsReal(DP)<AbsReal(DX) then
            begin
                INFO := 3;
                BOUND := True;
                THETA := 3*(FX-FP)/(STP-STX)+DX+DP;
                S := Max(AbsReal(THETA), Max(AbsReal(DX), AbsReal(DP)));
                GAMMA := S*SQRT(Max(0, Sqr(THETA/S)-DX/S*(DP/S)));
                if STP>STX then
                begin
                    GAMMA := -GAMMA;
                end;
                P := GAMMA-DP+THETA;
                Q := GAMMA+(DX-DP)+GAMMA;
                R := P/Q;
                if (R<0) and (GAMMA<>0) then
                begin
                    STPC := STP+R*(STX-STP);
                end
                else
                begin
                    if STP>STX then
                    begin
                        STPC := STMAX;
                    end
                    else
                    begin
                        STPC := STMIN;
                    end;
                end;
                STPQ := STP+DP/(DP-DX)*(STX-STP);
                if BRACKT then
                begin
                    if AbsReal(STP-STPC)<AbsReal(STP-STPQ) then
                    begin
                        STPF := STPC;
                    end
                    else
                    begin
                        STPF := STPQ;
                    end;
                end
                else
                begin
                    if AbsReal(STP-STPC)>AbsReal(STP-STPQ) then
                    begin
                        STPF := STPC;
                    end
                    else
                    begin
                        STPF := STPQ;
                    end;
                end;
            end
            else
            begin
                INFO := 4;
                BOUND := False;
                if BRACKT then
                begin
                    THETA := 3*(FP-FY)/(STY-STP)+DY+DP;
                    S := Max(ABSReal(THETA), Max(ABSReal(DY), ABSReal(DP)));
                    GAMMA := S*SQRT(Sqr(THETA/S)-DY/S*(DP/S));
                    if STP>STY then
                    begin
                        GAMMA := -GAMMA;
                    end;
                    P := GAMMA-DP+THETA;
                    Q := GAMMA-DP+GAMMA+DY;
                    R := P/Q;
                    STPC := STP+R*(STY-STP);
                    STPF := STPC;
                end
                else
                begin
                    if STP>STX then
                    begin
                        STPF := STMAX;
                    end
                    else
                    begin
                        STPF := STMIN;
                    end;
                end;
            end;
        end;
    end;
    if FP>FX then
    begin
        STY := STP;
        FY := FP;
        DY := DP;
    end
    else
    begin
        if SGND<0.0 then
        begin
            STY := STX;
            FY := FX;
            DY := DX;
        end;
        STX := STP;
        FX := FP;
        DX := DP;
    end;
    STPF := Min(STMAX, STPF);
    STPF := Max(STMIN, STPF);
    STP := STPF;
    if BRACKT and BOUND then
    begin
        if STY>STX then
        begin
            STP := Min(STX+0.66*(STY-STX), STP);
        end
        else
        begin
            STP := Max(STX+0.66*(STY-STX), STP);
        end;
    end;
end;


subroutine LBFGSNewIteration(const X : TReal1DArray;
     F : Double;
     const G : TReal1DArray);
begin
end;

end module lbfgs